#define c 137.037
#define pi 3.1416
#define nplim 2
#include <math.h>
#include <stdio.h>
#include <vector.h>
#include "../include/gridobject.h"
#include "../include/laserfield.h"

void CalculaRo(vector<double> &_xp, grid_object &_ro,vector<double> &_qp,double _deltag, int _na, FILE *fp0) //na is a parameter of the sheet
{
	//int np,npion,ndg,
	//int i,j,na,l;
	//double xp[np],Ro[ndg],xg,qp[np],deltag,deltat,Tmax,qm;
	double xa,xb,qa,qb,qi,qil;
	double xg;
	double qtotalfinal, qtotalinicial;
	
	qtotalinicial=0;
    for(int j=0;j<_xp.size();j++)
        qtotalinicial+=_qp[j];
    //printf("***Carga total antes de calculaRo= %.4f\n", qtotalinicial);
	
	for (int i=0;i<_ro.v.size();i++)
    {
		_ro.v[i]=0.;
    }
	
	for(int j=0;j<_xp.size();j++)
    {
		xg=_xp[j];
         
		//int i=int(xg/_deltag)+1;
        int i=int(xg/_deltag); //OJO!! MODIFICACION RESPECTO A FORTRAN, LE QUITO UNA UNIDAD PARA HACER EL GRID IGUAL
		//fprintf(fp0, "%d, %e, %e, %e\n", i, xg, xg/_deltag, (xg/_deltag)+1);
		//printf("\nPosicion=%.3f\tPunto en la red (i)=%d\tCarga=%.4f\n", xg, i, _qp[j]);
		//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
        qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
		qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
		//printf("qa=%.4f\tb=%.4f\n", qa, qb);
		
		xa=(xg+(i)*_deltag)/2.;
		xb=(xg+(i+1)*_deltag)/2.;
		qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
		
		for (int l=i-_na/2;l<=i+_na/2;l++)
		{
			if(l>=1 && (l<=_ro.v.size()))
			{	_ro.v[l]=_ro.v[l]+qi/_na/_deltag;
                //printf("qi:Asigno %.4f al punto %d\n", qi/_na/_deltag, l);
            }
		}
		qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
		//printf("qi=%.4f\tqi1=%.4f\n", qi, qil);
		for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
		{
			if( (l>=1) && (l<=_ro.v.size()))
			{	_ro.v[l]=_ro.v[l]+qil/_na/_deltag;
				//printf("qil:Asigno %.4f al punto %d\n", qil/_na/_deltag, l);
            }
		}
		
	}
	

    /*for (int i=0;i<_ndg;i++)
        fprintf(fp0, "[%d]%e\t",i,_ro[i]);
    fprintf(fp0, "\n");*/
    
    /*printf("***Carga total tras calculaRo= %.4f\n", qtotalfinal*_deltag);*/
} // Calculate Ro


void CalculaJ(int _label, vector<double> &_xp,vector<double> &_vpx, vector<double> &_vpy,vector<double> &_vpz,grid_object &_Jx,grid_object &_Jy,grid_object & _Jz,vector<double> &_qp,double _deltat,double _deltag, int _na, FILE *fp0, FILE *fp1)
{
//	int np,npion,ndg,i,j,na.l;
//	double xp[np],vpx[np],vpy[np],vpz[np];
//	double Jx[ndg],Jy[ndg],Jz[ndg],xg,qp[ndg],deltag,deltat,Tmax,gamma,qm;
	double xa,xb,qa,qb,qi,qil,q,mass;
	double xg;
	
	if(_label==1)
	{
		for(int i=0;i<_Jx.v.size();i++)
			_Jx.v[i]=0.;
		
		for(int j=0;j<_xp.size();j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/lightC_au/lightC_au);
			
            //MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_Jx.v.size()))
					_Jx.v[l]=_Jx.v[l]+qi/_na/_deltag*(_vpx[j])/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_Jx.v.size()))
					_Jx.v[l]=_Jx.v[l]+qil/_na/_deltag*(_vpx[j])/gamma;
			}
			//fprintf(fp1, "%d\t%e\t%e\t%e\n",j, _vpx[j], _vpy[j], _vpz[j]);
			
		}//end j
		
    /*for (int i=0;i<_ndg;i++)
        fprintf(fp0, "[%d]: %e\n",i,_Jx[i]);*/
        
    //fprintf(fp0, "\n");		
       
		
	}//End lable1=Jx
	
	
	
	
	if(_label==2)
	{
		for(int i=0;i<_Jy.v.size();i++)
			_Jy.v[i]=0.;
		
		for(int j=0;j<_xp.size();j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/lightC_au/lightC_au);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_Jy.v.size()))
					_Jy.v[l]=_Jy.v[l]+qi/_na/_deltag*_vpy[j]/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_Jy.v.size()))
					_Jy.v[l]=_Jy.v[l]+qil/_na/_deltag*(_vpy[j])/gamma;
			}
			
		}//end j
		
	}//end label2=Jy
	
	
	if(_label==3)
	{
		for(int i=0;i<_Jz.v.size();i++)
			_Jz.v[i]=0.;
		
		for(int j=0;j<_xp.size();j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/lightC_au/lightC_au);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=0) && (l<=_Jz.v.size()))
					_Jz.v[l]=_Jz.v[l]+qi/_na/_deltag*(_vpz[j])/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=0) && (l<=_Jz.v.size()))
					_Jz.v[l]=_Jz.v[l]+qil/_na/_deltag*(_vpz[j])/gamma;
			}
			
		}//end j
		
	}//end label==3
	
}//End CalculaJ


double E0x(double _xg, double _t)
{
	//E0x=0.01*sin(2.*w0*t)*sin(wo*t)*sin(w0*t);
	double ex=0.;
	return ex;
}

double E0y(double _xg,double _t,double _xgleft, laser & _laser0)
{
	double xx=_xgleft;
	double ey;
	
	if((_xg-lightC_au*_t < xx) && (_laser0.kp*(_xg-xx)-_laser0.wp*_t >=-Pi) )
		ey=_laser0.E0*sin(_laser0.kp*(_xg-xx)-_laser0.wp*_t)*sin(_laser0.kp*(_xg-xx)-_laser0.wp*_t)*sin(_laser0.k0*(_xg-xx)-_laser0.w0*_t);
	else
		ey=0.;
	return ey;
}

double E0z(double _xg, double _t)
{
	double ez=0.;
	return ez;
}

double B0x(double _xg,double _t)
{
	double bx=0.;
	return bx;
}

double B0y(double _xg,double _t)
{
	double by=0.;
	return by;
}

double B0z(double _xg,double _t,double _xgleft, laser _laser0)
{
	double bz=E0y(_xg, _t, _xgleft, _laser0);
	return bz;
}


void CalculaEx(grid_object & _Ex, grid_object & _Phi, grid_object & _DAx, double _t,double _deltat, FILE *fp0)
{
	// -grad Phi -1/c dAx/dt
	int ndg=_Phi.v.size();
	double deltag=_Phi.dx1;
	
	for(int i=1;i<(ndg-1);i++)
	{
		//_Ex[i]=-(_Phi[i+1]-_Phi[i-1])/2./_deltag-_DAx[i]/c/_deltat+E0x((i-1)*_deltag,_t);//MODIFICACION DE FORTRAN A C
		_Ex.v[i]=-(_Phi.v[i+1]-_Phi.v[i-1])/2./deltag-_DAx.v[i]/lightC_au/_deltat+E0x((i)*deltag,_t);
	}
	_Ex.v[0]=-(_Phi.v[1]-_Phi.v[0])/deltag-_DAx.v[0]/lightC_au/_deltat+E0x(0.,_t);
	_Ex.v[ndg-1]=-(_Phi.v[ndg-1]-_Phi.v[ndg-2])/deltag-_DAx.v[ndg-1]/lightC_au/_deltat+E0x((ndg-1)*deltag,_t);
	
	/*for(int j=0; j<_ndg; j++)
	         fprintf(fp0, "%d\t%e\n", j, _Ex[j]);*/
}


void CalculaEy(grid_object & _Ey,grid_object & _DAy,double _t,double _deltat, double _xg,double _xgleft, laser &_laser0)// double _wp,double _kp,double _E0,double _k0,double _w0)
{
	// -1/c dAx/dt
	int ndg=_DAy.v.size();
	double deltag=_DAy.dx1;
	
	for(int i=0;i<ndg;i++)
	{
		//_Ey[i]=-_DAy[i]/c/_deltat+ E0y((i-1)*_deltag,_t,_xgleft, _wp,_kp, _E0,_k0, _w0); MODIFICACION PARA C
		_Ey.v[i]=-_DAy.v[i]/lightC_au/_deltat+ E0y((i)*deltag,_t,_xgleft,_laser0);// _wp,_kp, _E0,_k0, _w0);
		//E0y( (i-1)*_deltag ,_t );
	}
}

void CalculaEz(grid_object & _Ez,grid_object & _DAz, double _t,double _deltat)
{
	// -1/c dAz/dt
	int ndg=_DAz.v.size();
	double deltag=_DAz.dx1;
	for(int i=0;i<ndg;i++)
	{
		//_Ez[i]=-_DAz[i]/c/_deltat+E0z((i-1)*_deltag,_t); //MODIFICACION PARA C
		_Ez.v[i]=-_DAz.v[i]/lightC_au/_deltat+E0z((i)*deltag,_t);
	}
}

void CalculaBy(grid_object & _By,grid_object & _Az, double _t,double _deltat)
{
	//-dAz/dx
	int ndg=_Az.v.size();
	double deltag=_Az.dx1;
	for(int i=1;i<(ndg-1);i++)
	{
		//_By[i]=-(_Az[i+1]-_Az[i-1])/2./_deltag+B0y((i-1)*_deltag,_t); MODIFICACION PARA C
		_By.v[i]=-(_Az.v[i+1]-_Az.v[i-1])/2./deltag+B0y((i)*deltag,_t);
	}
	_By.v[0]=-(_Az.v[1]-_Az.v[0])/deltag+B0y(0.,_t);
	//_By[_ndg-1]=-(_Az[_ndg-1]-_Az[_ndg-2])/_deltag+B0y((_ndg-1)*_deltag,_t); //MODIFICADO DE FORTRAN A C lo puse a 2, pero lo vuelvo a 1
	_By.v[ndg-1]=-(_Az.v[ndg-1]-_Az.v[ndg-2])/deltag+B0y((ndg-1)*deltag,_t); 
}

void CalculaBz(grid_object & _Bz,grid_object & _Ay,double _t,double _deltat,double _xg,double _xgleft, laser & _laser0)//double _wp,double _kp,double _E0,double _k0,double _w0)
{
	//dAy/dx
	int ndg=_Ay.v.size();
	double deltag=_Ay.dx1;
	for(int i=1;i<(ndg-1);i++)
	{
		//_Bz[i]=-(_Ay[i+1]-_Ay[i-1])/2./_deltag+B0z((i-1)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0); MODIFICACION PARA C
        _Bz.v[i]=(_Ay.v[i+1]-_Ay.v[i-1])/2./deltag+B0z((i)*deltag, _t, _xgleft, _laser0);// _wp, _kp, _E0, _k0, _w0);
   	}
	_Bz.v[0]=(_Ay.v[1]-_Ay.v[0])/deltag+B0z(0.,_t, _xgleft,_laser0);
	//_Bz[_ndg-1]=-(_Ay[_ndg-1]-_Ay[_ndg-2])/_deltag+B0z((_ndg-1)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0);//((_ndg-1),_t); MODIFICADO DE FORTRAN A C lo puse a 2 pero lo vuelvo a 1
	_Bz.v[ndg-1]=(_Ay.v[ndg-1]-_Ay.v[ndg-2])/deltag+B0z((ndg-1)*deltag, _t, _xgleft, _laser0);// _wp, _kp, _E0, _k0, _w0);
	
}


void interpola(double _x,grid_object & _Ex,grid_object & _Ey,grid_object & _Ez,grid_object & _By,grid_object & _Bz,double &_Epx,double &_Epy,double &_Epz,double &_Bpx,double &_Bpy,double &_Bpz, double _t,int _ndg,double _deltag, int _na, FILE *fp0)
{
	double xx=_x;
	//int i=int(xx/_deltag)+1;//MODIFICACION DE C
	int i=(int)round(xx/_deltag);
	_Epx=0.;
	_Epy=0.;
	_Epz=0.;
	_Bpx=0.;
	_Bpy=0.;
	_Bpz=0.;
	//if(i<309) printf("xx=%f, deltag=%f, xx/deltag=%f, (xx/deltag)+1=%f, int(xx/deltag)+1=%d\n", xx, _deltag, xx/_deltag, (xx/_deltag)+1, (int)(xx/_deltag)+1);
	//fprintf(fp0, "%d\t%e\txx/deltag=%e\t%e\t%e\t%e\t%e\t%e\t",i,_x, _x/_deltag, _Ex[i], _Ey[i], _Ez[i], _By[i], _Bz[i]);
	
	//printf("Int:");
	for(int l=i-_na/2;l<=i+_na/2;l++)
	{
		if((l>=0) && (l<_ndg) )
		{
			_Epx=_Epx+_Ex.v[l]/_na;
			_Epy=_Epy+_Ey.v[l]/_na;
			_Epz=_Epz+_Ez.v[l]/_na;
			
			_Bpx=_Bpx+B0x(xx,_t)/_na;
			_Bpy=_Bpy+_By.v[l]/_na;
			_Bpz=_Bpz+_Bz.v[l]/_na;
		}
	}
	//fprintf(fp0, "%e\t%e\t%e\t%e\t%e\t%e\n",_Epx, _Epy, _Epz, _Bpx, _Bpy, _Bpz);
	//printf("Epx=%.3f\tEpy=%.3f\tEpz=%.3f\tBpx=%.3f\tBpy=%.3f\tBpz=%.3f\t", *_Epx, *_Epy, *_Epz, *_Bpx, *_Bpy, *_Bpz);

}//end interpola

void CalculaPotenciales(grid_object & _phi, grid_object & _Ax, grid_object & _Ay, grid_object & _Az, grid_object & _Axold, grid_object & _Ayold, grid_object & _Azold, grid_object & _DAx, grid_object & _DAy, grid_object & _DAz, grid_object & _Jx, grid_object & _Jy, grid_object & _Jz,grid_object & _IntJx, grid_object & _IntJy, grid_object & _IntJz, grid_object & _Ro, grid_object & _IntRo, grid_object & _Rop, grid_object & _Roold, grid_object & _Ropold, grid_object & _IAxM, grid_object & _IAyM, grid_object & _IAzM, grid_object & _IAxP, grid_object & _IAyP, grid_object & _IAzP, grid_object & _IPhiM, grid_object & _IPhiP, double _deltat, FILE *fp0)
{
        int ndg=_Jx.v.size();
        double deltag=_Jx.dx1;
        
        for(int i=1;i<ndg-1;i++)
		{
			//Time integrals of Ro and J
			
			_IntRo.v[i]=_IntRo.v[i]+(_Ro.v[i]+_Rop.v[i]+_Roold.v[i]+_Ropold.v[i])*0.5*_deltat;  //Ros and j are defined in diff time positions
			
			_IntJx.v[i]=_IntJx.v[i]+_Jx.v[i]*_deltat;
			_IntJy.v[i]=_IntJy.v[i]+_Jy.v[i]*_deltat;
			_IntJz.v[i]=_IntJz.v[i]+_Jz.v[i]*_deltat;
			
			_phi.v[i]=dospi*lightC_au*(_IPhiM.v[i-1]+_IPhiP.v[i+1]+deltag*_IntRo.v[i]);
			
			_Ax.v[i]=dospi*(_IAxM.v[i-1]+_IAxP.v[i+1]+deltag*_IntJx.v[i] );
			_Ay.v[i]=dospi*(_IAyM.v[i-1]+_IAyP.v[i+1]+deltag*_IntJy.v[i] );
			_Az.v[i]=dospi*(_IAzM.v[i-1]+_IAzP.v[i+1]+deltag*_IntJz.v[i] );
		}
		
		_IntRo.v[0]=_IntRo.v[0]+(_Ro.v[0]+_Rop.v[0]+_Roold.v[0]+_Ropold.v[0])*0.5*_deltat;
		
		_IntJx.v[0]=_IntJx.v[0]+_Jx.v[0]*_deltat;
		_IntJy.v[0]=_IntJy.v[0]+_Jy.v[0]*_deltat;
		_IntJz.v[0]=_IntJz.v[0]+_Jz.v[0]*_deltat;
		
		_phi.v[0]=dospi*lightC_au*(_IPhiP.v[1]+deltag*_IntRo.v[0]);
		_Ax.v[0]=dospi*(_IAxP.v[1]+deltag*_IntJx.v[0]);
		_Ay.v[0]=dospi*(_IAyP.v[1]+deltag*_IntJy.v[0]);
		_Az.v[0]=dospi*(_IAzP.v[1]+deltag*_IntJz.v[0]);
		
		_IntRo.v[ndg-1]=_IntRo.v[ndg-1]+(_Ro.v[ndg-1]+_Rop.v[ndg-1]+_Roold.v[ndg-1]+_Ropold.v[ndg-1])*0.5*_deltat;
        _IntJx.v[ndg-1]=_IntJx.v[ndg-1]+_Jx.v[ndg-1]*_deltat;
		_IntJy.v[ndg-1]=_IntJy.v[ndg-1]+_Jy.v[ndg-1]*_deltat;
		_IntJz.v[ndg-1]=_IntJz.v[ndg-1]+_Jz.v[ndg-1]*_deltat;
		
		_phi.v[ndg-1]=dospi*lightC_au*(_IPhiP.v[ndg-2]+deltag*_IntRo.v[ndg-1]);
		_Ax.v[ndg-1]=dospi*(_IAxP.v[ndg-2]+deltag*_IntJx.v[ndg-1]);
		_Ay.v[ndg-1]=dospi*(_IAyP.v[ndg-2]+deltag*_IntJy.v[ndg-1]);
		_Az.v[ndg-1]=dospi*(_IAzP.v[ndg-2]+deltag*_IntJz.v[ndg-1]);
		
		for (int i=0;i<ndg;i++)
		{
			_DAx.v[i]=(_Ax.v[i]-_Axold.v[i]);
			_DAy.v[i]=(_Ay.v[i]-_Ayold.v[i]);
			_DAz.v[i]=(_Az.v[i]-_Azold.v[i]);			
		}
		
		//Spatial integrals.
		for (int i=0;i<ndg-1;i++)
		{
			_IPhiP.v[i]=_IPhiP.v[i+1]+_IntRo.v[i]*deltag;
			_IAxP.v[i]=_IAxP.v[i+1]+_IntJx.v[i]*deltag;
			_IAyP.v[i]=_IAyP.v[i+1]+_IntJy.v[i]*deltag;
			_IAzP.v[i]=_IAzP.v[i+1]+_IntJz.v[i]*deltag;
		}
		
		_IPhiP.v[ndg-1]=_IntRo.v[ndg-1]*deltag;
		_IAxP.v[ndg-1]=_IntJx.v[ndg-1]*deltag;
		_IAyP.v[ndg-1]=_IntJy.v[ndg-1]*deltag;
		_IAzP.v[ndg-1]=_IntJz.v[ndg-1]*deltag;
		
		for(int i=ndg-1;i>=1;i--)
		{
			_IPhiM.v[i]=_IPhiM.v[i-1]+_IntRo.v[i]*deltag;
			_IAxM.v[i]=_IAxM.v[i-1]+_IntJx.v[i]*deltag;
			_IAyM.v[i]=_IAyM.v[i-1]+_IntJy.v[i]*deltag;
			_IAzM.v[i]=_IAzM.v[i-1]+_IntJz.v[i]*deltag;
		}

		_IPhiM.v[0]=_IntRo.v[0]*deltag;
		_IAxM.v[0]=_IntJx.v[0]*deltag;
		_IAyM.v[0]=_IntJy.v[0]*deltag;
		_IAzM.v[0]=_IntJz.v[0]*deltag;
		
		for(int i=0;i<ndg;i++)
		{
			_Axold.v[i]=_Ax.v[i];
			_Ayold.v[i]=_Ay.v[i];
			_Azold.v[i]=_Az.v[i];

			_Roold.v[i]=_Ro.v[i];
			_Ropold.v[i]=_Rop.v[i];
		}
}






void ParticleMover(double &_xp, double &_vpx,double &_vpy, double &_vpz, double _Epx,double _Epy,double _Epz,double _Bpx,double _Bpy,double _Bpz, double _deltat, double ae)
{	
            double vvx,vvy,gvxs,gvys,vvz,gvzs,dcg,fb;
            double aexpt,aeypt,aezpt,abxpt,abypt,abzpt;
            double gamma;
    
           	vvx=0.;
        	vvy=0.;
        	vvz=0.;
        	
        	gvxs=0.;
        	gvys=0.;
        	gvzs=0.;
        	
        	aexpt=0.;
        	aeypt=0.;
        	aezpt=0.;
    
            gamma=sqrt(1.+(_vpx*_vpx+_vpy*_vpy+_vpz*_vpz)/lightC_au/lightC_au );
			
			
			aexpt=ae*_Epx;
			aeypt=ae*_Epy;
			aezpt=ae*_Epz;
			
			abxpt=ae*_Bpx;
			abypt=ae*_Bpy;
			abzpt=ae*_Bpz;
			
			gvxs=_vpx+aexpt;
			gvys=_vpy+aeypt;
			gvzs=_vpz+aezpt;
			
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
			dcg=1./c/gamma;
			
			abxpt=dcg*abxpt;
			abypt=dcg*abypt;
			abzpt=dcg*abzpt;
			
			vvx=gvxs+gvys*abzpt-gvzs*abypt;
			vvy=gvys+gvzs*abxpt-gvxs*abzpt;
			vvz=gvzs+gvxs*abypt-gvys*abxpt;
			//fprintf(fp11, "%e %e %e ", vvx, vvy, vvz);
			//Second half rotation and second half acceleration.
			fb=2./(1.+abxpt*abxpt+abypt*abypt+abzpt*abzpt);
			abxpt=abxpt*fb;
			abypt=abypt*fb;
			abzpt=abzpt*fb;
			
			gvxs=gvxs+vvy*abzpt-vvz*abypt+aexpt;
			gvys=gvys+vvz*abxpt-vvx*abzpt+aeypt;
			gvzs=gvzs+vvx*abypt-vvy*abxpt+aezpt;
			gamma=sqrt(1.+(gvxs*gvxs+gvys*gvys+gvzs*gvzs)/lightC_au/lightC_au);
			
			_vpx=gvxs;
			_vpy=gvys;
			_vpz=gvzs;
			
			_xp=_xp+_deltat*_vpx/gamma;

}

void CalculaCamposInternos(grid_object & _Ex, grid_object & _Ey, grid_object & _Ez, grid_object & _By, grid_object & _Bz, grid_object & _Jx, grid_object & _Jy, grid_object & _Jz, grid_object & _Fp, grid_object & _Fm, grid_object & _Gp, grid_object & _Gm,  double _t,double _deltat, double _xgleft, laser _laser0, FILE *fp0)
{
     
    int ndg=_Ex.v.size();
    double deltag=_Ex.dx1;
    double fp, fp_old, gm, gm_old, fm, fm_old, gp, gp_old;
    
    fp_old=_Fp.v[0];
    gm_old=_Gm.v[0];
    
    //Condiciones para la primera posicion del grid
    
    /***********************************************/
    //_Fp.v[0]= (E0y(0,_t,_xgleft, _laser0)+B0z(0,_t,_xgleft, _laser0))/2.;
    //_Gm.v[0]= (E0z(0,_t)+B0y(0,_t))/2.;
    
    _Fp.v[0]= 0;
    _Gm.v[0]= 0;    
    
    //cell->fp = laser_front.Qy * laser_front.field( time );
    //cell->gm = laser_front.Qz * laser_front.field( time + laser_front.shift );
    /***********************************************/                                   
    
    _Fm.v[0]=_Fm.v[1]- 2.*Pi*_deltat*_Jy.v[1];
    _Gp.v[0]=_Gp.v[1]- 2.*Pi*_deltat*_Jz.v[1];
    
    _Ey.v[0]=_Fp.v[0]+_Fm.v[0];
    _Bz.v[0]=_Fp.v[0]-_Fm.v[0];
    _Ez.v[0]=_Gp.v[0]+_Gm.v[0];
    _By.v[0]=_Gp.v[0]-_Gm.v[0];
    
    _Ex.v[0]= _Ex.v[0]-4.*Pi*_deltat*_Jx.v[0];
          
          
    for(int i=1; i<(ndg-1); i++)
    {
          fp=_Fp.v[i];
          _Fp.v[i]= fp_old- 2.*Pi*_deltat * _Jy.v[i-1];
          fp_old=fp;
          
          gm=_Gm.v[i];
          _Gm.v[i]= gm_old- 2.*Pi*_deltat*_Jz.v[i-1];   
          gm_old=gm;                                       
          
          _Fm.v[i]=_Fm.v[i+1]- 2.*Pi*_deltat*_Jy.v[i];
          _Gp.v[i]=_Gp.v[i+1]- 2.*Pi*_deltat*_Jz.v[i];
    
          _Ey.v[i]=_Fp.v[i]+_Fm.v[i];
          _Bz.v[i]=_Fp.v[i]-_Fm.v[i];
          _Ez.v[i]=_Gp.v[i]+_Gm.v[i];
          _By.v[i]=_Gp.v[i]-_Gm.v[i];
          
          _Ex.v[i]= _Ex.v[i]-4.*Pi*_deltat*_Jx.v[i];
    }            
    
    
    //Condiciones para la ultima posicion del grid
    _Fp.v[ndg-1]= fp_old- 2.*Pi*_deltat*_Jy.v[ndg-2];
    _Gm.v[ndg-1]= gm_old- 2.*Pi*_deltat*_Jz.v[ndg-2];   
                                       
    /***********************************************/
    /*_Fm.v[ndg-1]= (E0y(ndg*_deltag,t,_xgleft, _laser0)-B0z(ndg*_deltag,t,_xgleft, _laser0))/2.;
    _Gp.v[ndg-1]= (E0z(ndg*_deltag,t)-B0y(ndg*_deltag,t))/2.;*/
    
    _Fm.v[ndg-1]= 0;
    _Gp.v[ndg-1]= 0;
    /***********************************************/
    
    _Ey.v[ndg-1]=_Fp.v[ndg-1]+_Fm.v[ndg-1];
    _Bz.v[ndg-1]=_Fp.v[ndg-1]-_Fm.v[ndg-1];
    _Ez.v[ndg-1]=_Gp.v[ndg-1]+_Gm.v[ndg-1];
    _By.v[ndg-1]=_Gp.v[ndg-1]-_Gm.v[ndg-1];
    
    _Ex.v[ndg-1]= _Ex.v[ndg-1]-4.*Pi*_deltat*_Jx.v[ndg-1];
    
}
/*
void Ionizacion(int _ndg, int _nplim,double _deltap0[],double _p0[],double _Ex[],double _Ey[],double _Ez[],double _wrate,double _t,double _deltat,int _ileft,int _iright,double _coef,double _exponente)//,int _ppp[][]
{
	double wrate[_ndg];
	for(int i=_ileft;i<=_iright;i++)
	{
		double aux=sqrt(_Ex[i]*_Ex[i]+_Ey[i]*_Ey[i]+_Ez[i]*_Ez[i]);
		if(aux>3.333e-3)
		{
			wrate[i]=_coef/aux*exp(-_exponente/aux);
		}
		else
			wrate[i]=0.;
		
		_deltap0[i]=_deltat*wrate[i]*_p0[i];
		if(_deltap0[i]>_p0[i] )
			_deltap0[i]=_p0[i];
		//_ppp[i][1]=0;
	}
}


void CalculaP0(double _xp[],double _vpx[],double _vpy[],double _vpz[],int _np,int _npion, int _ndg,double _qp[],int _nplim, double _ppp[][nplim+1],double _p0[],double _rop[],double _deltap0[], int n, int _iseed,double _deltag,int _ileft, int _iright,double _cargaminima,int _nppp,int _npminimo,int _npsinp)
{
	int jal=1;//+int(_npion*rand(_iseed));
	for(int j=jal;j<_npion;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j])< (_xp[j]-(i-1)*_deltag )    )
			i=i+1;
		if(_ppp[i][1] < _nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;//+_ppp[i][1]]=j;  //Check this line
		}
	}
	for(int j=1;j<jal-1;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j]) < (_xp[j]-(i-1)*_deltag)  )
		{
			i=i+1;
		}
		if(_ppp[i][1]<_nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;// +_ppp[i,1]]=j;  fix this !!!
		}
	}
	
	for(int i=_ileft;i<_iright;i++)
	{
		if(_deltap0[i]> _cargaminima)
			_nppp=_ppp[i][1];
		if(_nppp!=0 && _npion>_npminimo)
		{
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
			double aux=_deltap0[i]*_deltag/_nppp;
			for(int j=1;j<_nppp;j++)
			{
				int jb=_ppp[i][1+j];
				double aux2=_qp[jb]-aux;
				if(aux2!=0)
				{
					_vpx[jb]=_vpx[jb]*_qp[jb]/aux2;
					_vpy[jb]=_vpy[jb]*_qp[jb]/aux2;
					_vpz[jb]=_vpz[jb]*_qp[jb]/aux2;
					_qp[jb]=aux2;
				}
			}//end for
		}//end if
		else if(_npion > _np)
		{
			_npion=_npion+1;
			_qp[_npion]=-_deltap0[i]*_deltag;
			_vpx[_npion]=0.;
			_vpy[_npion]=0.;
			_vpx[_npion]=0.;
			_xp[_npion]=_deltag*(i-1.);
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
		}
		else
			_npsinp=_npsinp+1;
	}
	
}


double ran1(int _idum,int _iff)
{
	int m=714025;
	int ia=1366;
	int ic=150889;
	double rm=1./m;
	double temp;
	int ir[97];
	int iy;
	
	if(_idum > 0 || _iff==0)
	{
		_iff=1;
		//_idum=mod(ic-_idum,m);
		for(int j=1;j<97;j++)
		{
			//_idum=mod(ia*_idum+ic,m);
			ir[j]=_idum;
		}
		//_idum=imod[ia*_idum+ic,m];
		iy=_idum;
	}
	int j=1+(97*iy)/m;
	//if(j>97 || j<1)
	//	break;
	
	iy=ir[j];
	double tempran1=iy*rm;
	//temp=imod(ia*_idum+ic,m);
	ir[j]=_idum;
	return tempran1;
}
*/

/*
void realft(double data,int n,int isignn)
{
	double pi=4.*atan(1.);
	double theta=pi/n;
	double cl=0.5;
	
	if(isign==1)
	{
		c2=-0.5;
		four1(data,n,1);
	}
	else
	{
		c2=0.5;
		theta=-theta;
	}
	double wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
	double wpi=sin(theta);
	double wr=1.+wpr;
	double wi=wpi;
	n2p3=2*n+3;
	for(int i=2;i<n/2;i++)
	{
		int il=2*i-1;
		int i2=il+1;
		int i3=n2p3-i2;
		int i4=i3+1;
		double wrs=wr;
		double wis=wi;
		
		double hlr=c1*(data[i1]+data[i3]);
		double hli=cl*(data[i2]-data[i4]);
		double h2r=-c2*(data[i2]+data[i4]);
		double h2i=c2*(data[i1]-data[i3]);
		
		data[i1]=h1r+wrs*h2r-wis*h2i;
		data[i2]=h1i+wrs*h2i+wis*h2r;
		data[i3]=h1r-wrs*h2r+wis*h2i;
		data[i4]=-h1i+wrs*h2i+wis*h2r;
		double wtemp=wr;
		
		double wr=wr*wpr-wi*wpi*wr;
		double wi=wi*wpr+wtemp*wpi+wi;
	}
	
	if(isign==1)
	{
		h1r=data[1];
		data[1]=h1r+data[2];
		data[2]=h1r-data[2];
	}
	else
	{
		h1r=data[1];
		data[1]=c1*(h1r+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n,-1);
	}
}


void four1(double data,int nn, int isign)
{
	double pi=4.*atan(1.);
	int n=2*nn;
	int j=1;
	
	for(int i=0;i<n;i++)
	{
			if(j>i)
			{
				tempr=data[j];
				tempi=data[j+1];
				
				data[j]=data[i];
				data[j+1]=data[i+1];
				data[i]=tempr;
				data[i+1]=tempi;
			}
			m=n/2;
			if(m>=2 && j > m)
			{
				j=j-m;
				m=m/2;
			}
			j=j+m;
			mmax=2;
			if(n>mmax)
			{
				istep=2*mmax;
				theta=2.*pi/(isign*mmax);
				wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
				wpi=sin(theta);
				wr=1.;
				wi=0.;
				for(int m=0;m<mmax;m++2)
				{
					for(int i=m;i<n;i++step)
					{
						j=i+mmax;
						tempr=wr*data[j]-wi*data[j+1];
						tempi=wr*data[j+1]-wi*data[j];
						
						data[j]=data[i]-tempr;
						data[j+1]=data[i+1]-tempi;
						data[i]=data[i]+tempr;
						data[i+1]=data[i+1]+tempi;
					}
					wtemp=wr;
					wr=wr*wprwi*wpi+wr;
					wi=wi*wpr+wtemp*wpi+wi;
				}
				mmax=istep;
			}
	}
}

*/
